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Abstract: The finite volume element method is combined with the characteristic method to treat the 
one-dimensional semiconductor device simulation problem. A fully discrete implicit Euler 
characteristic finite volume element scheme is derived and studied. With different time steps 
for the electrostatic potential and other unknown quantities, the computational procedure 
of the new method is constructed. A priori convergence estimate of first-order accuracy in 
L?-norm is obtained under suitable assumptions. Numerical experiments demonstrate the 
validity of the theoretical results and efficiency of the method. 
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1 Introduction 


We consider the following one-dimensional semiconductor device simulation in Q = fa, b]: 


Ow 7 
-3 = (p -e+ N(2)), (x,t) € Q x (0,T], (1) 
ð o ð 0 = 
= = z (Del@) 5 T peee) - R(e,p, T), (æ,t) €Q x (0,7), (2) 
0; oO re) ð = 
a z zz (Pole) Dz + mp) - R(e,p,T), (x,t) € Q x (0,T], (3) 


ôT &@T 0 ð 
oa) — Sax = (oE + mla) pE) 


= (D.(@) & 5 pe(a)eo*) a (z,t)EQ2x (0, f}, (4) 
where (1) is the electrostatic potential equation, (2) and (3) are the electron and hole concen- 
tration equations, respectively, and (4) is the temperature equation. The four equations with 
initial and boundary conditions make up a closed system. œ = q/e, where q and £ are both 
positive constants, respectively, representing the electron charge and the dielectric permittiv- 
ity. D,(x)(s = e,p) is the diffusion coefficient and s(x) (s = e,p) is the mobility. D,(z) 
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and ps(z) obey the Einstein relation D,(z) = Urp.(x), where Ur is the thermal voltage. 
N(x) = Np(x) — Na(z) is a given function, where Np(z) is the donor impurity concentration 
and Na(x) is the acceptor impurity concentration. R(e,p,T) is a recombination term of the 
electron and hole concentrations with the influence of temperature. For (1)-(4), we consider 
the initial value conditions 


e(x,0) = eo(x), p(x,0)=po(x), T(z,0)=To(z), rE, (5) 


and the initial value of y can be calculated by (1) and (5). For convenience, we assume that 
problems (1)-(5) are 9-periodic!*], For ordinary boundary value problems we can obtain the 
similar results via dealing with the technique extension or mirror reflection!). 

Since Gummell4! first proposed sequence iterative computational methods for this kind of 


problem, there have been many numerical contributions on such problems|3:5-41, 


Among the 
semiconductor mathematical models, the electron concentration equation (2) and the hole con- 
centration equation (3) are always convection-dominated. It means that the convection term 
is distinctly more important than the diffusion term. If we only use pure standard methods 
such as the finite element method, the finite difference method or the finite volume element 
method, it will generate numerical dispersion and nonphysical oscillations unless the mesh is 
excessively impractical. So Douglas and Russelll’-§] proposed the characteristic method and 
combined the method with the finite element method or the finite difference method. They an- 
alyzed the method under Q-periodic assumptive condition and solved many physical problems 
with satisfactory results. 

In this paper, with the modified characteristic method, we construct a fully discrete implicit 
Euler characteristic finite volume element method for the one-dimensional semiconductor device 
simulation. With different time steps for the electrostatic potential and other unknown quan- 
tities, we construct the computational procedure of the method. A priori convergence estimate 
of first-order accuracy in the L?-norm is obtained under some suitable assumptions. In the end 


we give numerical experiments to validate the theoretical results. 


2 Characteristic finite volume element scheme 


Let us discrete the domain 2 = [a,b] into J subintervals by the points a = rp < 41 < T2 < 
may <i Say =b. Let hi = ti — 24-1, I, = [zi-1, zi] and h = maxi<igy, 1=1,2,--+ J. 
Assume that the partition is quasi-uniform, i.e., there is a positive constant Co > 0 such that 
hi > Coh for all i = 1,2,+++ , J. This partition is denoted by Th = {J;}7_, and the subinterval 
I; is called finite element. 

We introduce the dual partition Tf with nodes 


Ti-1/2 = (ti-1+4;)/2, 1<i<J. 


Let 
I = [20,212], Tf = [@i-1/2, Zisiye], 13 = [@s-1/2, tg], LSt<J-1. 


I} is also called control volume. Let 


H+ (Q) = {u(x) € H! (9), Q — periodic}, U = A! (N) N H?(Q). 
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Then we introduce the solution space Ua C U of the piecewise linear functions on the original 


partition of the domain. The basis function with respect to zi isl] 


1> htr- zil, Ti-1 STX Ti, 





gi(z)=4 1- hghle-— zil, Ti LTL Tiy, 


0, otherwise. 


J 
Thus, for any up € Un, we have up = >> uiy;(x), where u; = up(zi). 
i= 


We choose the test function space Vp of piecewise constant functions on the dual partition 
T}. The characteristic function, defined by 


> l, TZi—1/2 Š T < Lire, 
pi (x) = A 
0, otherwise, 


J 
forms a basis for the space Va. For any up € Vh, there is vn = >> vipž (x), where v; = va (zi). 
i=l 
By means of the local basis we can introduce the interpolation operators: Ip : U — Up, and 


If :U — Vn such that 


J J 


nu = >X ulzi)el), Tu=) ulel), wel. 


i=1 t=1 


Define the bilinear form 


D—— dz, (u,v) EUn x Un, 


23y 8 
=, z (D5")vdr, (u,v) EUn X Va, 
and a(l; un, vh) = alun, vn). To give a fully discredited scheme, we take the time step At = 
T/N, where N is a positive integer and t? = nAt, n = 0,1,---, N. Denote u(t") = uf. 

Let Te(£) be the unit vector of (— Heu, 1) and 7p(z) be the unit vector of (upu, 1), respectively. 
Denote ¢, = [1 + y2u?]!/? and s = e,p. Then the characteristic derivative is approximated by 





ð -1(2- u) oat (Fy 2 ) 
OTe be \Ot He Dad” OT, bp \Ot Petar) 
We use a backward difference quotient for ge" (x) = 2 (z, t”) along the characteristics. Specif- 


ically, we take 


Bel eine Se (adie) E ea Lew DY) 
Or. ea Ae T Ai 





$e(z) 3 (6) 

Since the electric potential on time steps always varies slowly in the semiconductor model, 
we will use large time steps to compute. However, we use small time steps to compute the 
electron concentration, the hole concentration and the temperature equations. We introduce 
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the following notation: At—the concentration and temperature time step, At®°—the first time 
step of the electric potential equation, and At—the other time step of the electric potential 
equation. Let 


j=At/At, j? =At°/At, t?=nAt, tm = Af? + (m—1)At. 
For a function w(z, ¢), let 
w= P(t"), Ym =V(tm), Hb" = (bt? — p")/At, 
Opm = (Wm+i — bm)/At (m>0), Abo = (di — 4o)/At®, 
and 
po, t <t, 
Ey? =2 (1+ Da -— Two, tı <t” < ty, t” = tı + VAT, (7) 
(1+ +) em — 5 Pat tm < t” < tmy, t =tmt VAt. 


The subscript is corresponding to the electric potential time level and the superscript is corre- 
sponding to the concentration and temperature time level, respectively. Ey", which is made 
up of the previous two time level values, denotes the extrapolation of ~ at t”. Let 


bh get 


Uh = An” Le z + peEuRAt, en = en(ee*); 
os =x-—ppEupAt, pe = palz"). 


Then a fully-discrete implicit Euler characteristic FVE scheme for the 1-D semiconductor prob- 
lem can be defined as to find {Yh,m, ek, pe, Tg} € Un, for n = 1,2,--- , N, such that 


a(Ph,m: Xh) = O(Phym — enm + Nm, Xn), VXn © Vn; (8) 


=n—1 
ep —2 n =n— nô e 
(o ae wn) + a(De; eh, Wr) — (e; ‘Euk wn) 


— (apeeR (PR! — ER! +N”) wa) + (R(eR pR, TR) wh) =0, Vun EVn, (9) 


(7 -pp 
At 


+ (appph (DR — ER! + N”) wn) + (R( eR, pR, TR), wn) =0, Wwn © Vr, (10) 


$ ð 
,wn) + a( Dp; Dh wh) + (r Eup? wn) 


ye i ala Oph. n n n 
(o-2 =~. zn) +a(Tk, zh) = - (Dp Euk, zn) + (Uppy (EuR), zn) 
ð m 
+ (D-5 Buk, zn) + (ueer(EuR)?, zn), V zn E Va, (11) 
ep(z) = (Ineo)(x), palz) = (Inpo)(£), TR(z) = UnTo)(x), «EQ. (12) 


The algorithm of the computational procedure is as follows. The initial values of el, p}, T 


are defined by (12); from (8) we can get Yhn,o. Assuming that the solutions of {Yh,m-1, ef}; 
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pr-', TẸ} are obtained, then we can compute the values of ef, pp and Tẹ by (9), (10) and 
(11). Repeat j times in this way, we can obtain Wa,m from (8). The whole computational 
procedure of this algorithm is as follows. e®, p?, TP, Pho; eh, pL, TH, Ps e, pe, TF, WPh,1; 
P a piott TEH... efti, ploti, Ti, app a; +e 


3 


3 Error estimate analysis 


Now, we consider error estimates in the L?-norm for the fully discrete characteristic finite 
volume element scheme (8)-(12). We will introduce some elliptic projections. Let y = dn € Un, 
such that 


alp — P, xh) =0, Vxn € Vh, (13) 
and let é, p, T € Up, satisfy 
a(Deje — č wh) + Aele — E,wn) =0, Von € Vn, (14) 
a(Dp;p — P, wh) + Ap(p — B, wn) = 0, Vwn € Va, (15) 
a(T -T wn) +Ar(T ~T,wpn) =0, Won € Vn, (16) 


where \,(s = e, p, T) are positive constants. It is well known that for {y,e,p,T} € U, we have 
IIb — dlr < Chivlle, (17) 


-3 
ls — Sly + 2 -D < ch( llsil2 + |2 =): s=e,p, T. (18) 
From (8) and (13), we get the error estimate equation of the electric potential 
alm, Xh) = A(Ph,m — Pm + €m — Ehm: Xn), VXh E Va. (19) 
We choose xn = [fm as the test function in (19) and obtain 
lômlli < C{[l€emll + Ilr] + A}- (20) 


For e, p and T, with the help of the elliptic projections and some other estimate technique, 


under some suitable assumptions we can get 


Bn {lée l? + les? + Weel? bod (lee i+ lép li + IEZI At 


n=1 
< E E eer (21) 
By (17), (18), (20), (21) and the triangle inequality, we conclude the theorem. 
Theorem 1 Let {e, p, T, y} be the solution of (1)-(5) and {e, p, T, Y} EU, {e}, p}, TR, 
Wh,m} be the solution of the fully-discrete implicit Euler characteristic finite volume element 


scheme (8)-(12), respectively. When (h, At) is small enough and satisfies h = O(At), Af? = 
O((At)?/3) and At = O((At)!/?), we have a prior estimate in the L?-norm as 


sup {|le” — el] + lp" — pal + IT” — TRI} + sup [lm — Yamll 
1<n<N 1<m<K 


< C(h+ At + (At®)9/? + (Ad)?). (22) 


930 CHINESE JOURNAL OF ENGINEERING MATHEMATICS VOL. 27 





4 Numerical experiments 


We consider the following system of coupled partial differential equations 


o? ð 

-F = ae, rEN=(0,2], t>0, 
2 

O65 One OR aia eQ = [0,2], t>0 s 

ôt = Ox2 Ox ? , T = ’ ’ , 

with exact solutions 
y= et” e372? +452—18 u = —(—74r + 45)et e7372" +452—18 

(24) 


e(z, t) = | — 74 + (-74a + 45)?] et e7372 +452—18 


where f(z, t) is decided by (24). Let the space-step h = 0.025 and the time-step At = 0.05. We 
use the characteristic finite volume element scheme to solve system (23) and obtain the discrete 


solution Yp and e. The error estimates in the L?-norm of Y and e are shown in Table 1 and 
Table 2. 


Table 1: Error estimate in the L?-norm when «€ = 0.01 








L?-norm = 0.1 t=0.3 t=0.5 t= 0.7 
lY — Yall 0.000190 0.000527 0.001001 0.001814 
lle — en|| 0.002768 0.007080 0.010890 0.017622 





Table 2: Error estimate in the L?-norm when «e = 0.0001 











L?-norm t=0.1 t=0.3 t=0.5 t=0.7 
Ib — vall 0.000516 0.001184 0.002232 0.003646 
lle — enll 0.003120 0.008211 0.014636 0.026958 


From Table 1 and Table 2, we can see that the proposed method is useful and effective even 
when c€ is very small. 
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EEH Be RUER MEN Sp Ee Be BR AR TK 
R, Wate? 
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